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1  SUMMARY 

The  AFRL-VT-WSU-UMD  Collaborative  Center  in  Multidisciplinary  Sciences  (CCMS)  was 
established  after  a  successful  proposal  developed  jointly  by  Virginia  Tech  and  Wright  State 
University  and  another  one  from  the  University  of  Maryland. 

These  proposals  were  developed  in  response  to  a  Broad  Agency  Announcement  (BAA  PKV-08- 
09)  with  the  following  stated  purpose: 

The  Multidisciplinary  Technology  Center  (MDTC)  of  AFRL  will  set  up  an  ongoing  partnership 
by  establishing  a  Collaborative  Center  in  Multidisciplinary  Sciences  (CCMS).  The  best  fit  will  be 
a  research  team  whose  expertise  strongly  complements  the  AFRL  team,  as  this  collaborative 
effort  will  be  an  extension  of  AFRL’s  Multidisciplinary  Technology  Center.  Developing  a 
Collaborative  Center  is  expected  to  increase  the  agility  and  responsiveness  of  AFRL  research 
efforts. 

This  report  contains  work  carried  out  at  the  University  of  Maryland,  from  the  end  of  the  Fall 
2009  semester.  Our  focus  included  the  following: 

1 .  Develop  tools  for  the  study  and  design  of  micro-air-vehicles  (MAVs). 

2.  Analyze  the  physics  observed  through  experiments  to  determine  beneficial  phenomena. 

3.  Compare  the  utility  of  models  of  varying  levels  of  fidelity  to  see  how  they  affect  the  design 
process. 

In  Section  2,  the  academic  activities  of  the  students  are  listed,  including  a  list  of  publications. 
The  development  and  analysis  of  two-dimensional  (2D)  flapping  models  are  discussed  in  Section 
3.  Special  attention  is  given  to  comparison  of  results  from  the  DNS  and  UVLM  models.  Novel 
vibration  experiments  carried  out  with  the  wings  of  living  hawkmoths  are  included  in  Section  4, 
along  with  a  formulation  of  a  highly  flexible  three-dimensional  (3D)  solid  model  for  use  in  FSI. 
This  section  also  contains  a  discussion  of  the  fluid-structure  interaction  (FSI)  implementation  of 
the  FSI  algorithm  inside  high-performance  FLASH  code.  Numerical  demonstrations  of  the 
effectiveness  of  the  implementation  are  presented  in  Section  5.  Finally,  concluding  remarks  are 
collected  together  in  Section  6. 
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3  Two-Dimensional  Modeling 

As  an  initial  step  towards  understanding  phenomena  associated  with  flapping  wings  and  tools 
that  can  be  used  to  study  these  systems,  two-dimensional  (2D)  systems  are  considered.  In  a  2D 
setting,  the  complicated  three-dimensional  (3D)  motions  are  idealized  to  motions  of  a 
representative  cross-section  of  a  wing,  as  depicted  in  Figure  1.  This  simplification  greatly 
reduces  the  complexity  of  the  physics,  the  possible  parameter  space  for  exploration,  and  the 
associated  computational  costs.  By  working  with  2D  models,  the  investigators  gained  a  basis  to 
understand  the  types  of  models,  which  would  need  to  be  constructed  for  the  3D  cases  to  be 
studied  later. 

3.1  Initial  study  on  flexibility 

In  the  first  model  studied,  a  single  measure  for  the  chord-wise  deflection  is  used.  Here,  there  are 
two  rigid  elements  connected  at  point  b  by  a  torsion  spring,  as  depicted  in  Figure  2.  This  torsion 
spring  is  assumed  to  be  linear  with  respect  to  the  deflection  angle  a  between  the  two  rigid 
elements.  The  location  of  point  b  is  specified  by  the  coordinates  (x,y),  and  9  is  the  orientation 
of  link  B.  The  center  of  mass  ,  i  E  {A, B),  of  each  link  is  located  a  distance  rh  from  the 

connection  point  b.  The  length  of  the  profile  is  1. 


Figure  1:  Schematic  illustration  of  a  two-dimensional  wing  profile  as  a  representative  cross- 

section  of  an  insect  wing^ 


^  Photo  of  female  Villa  hottentotta,  used  under  a  Creative  Commons  Attribution- Share  Alike  license, 
http://commons.wikimedia.0rg/wiki/File:Villa  hottentotta  female. ipg 
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Figure  2:  Schematic  of  two-dimensional  profile  with  discrete  flexibility.  (A)  Details  of  the 
profile’s  description  and  coordinates  on  the  10%  thickness  case.  (B)  Example  of  the  3%  thickness 

case. 

The  equations  of  motion  for  this  structural  model  can  be  arranged  as 


mA+nig  0  sin  9  —  m/rjA  sin(Q!  +  6)  — 171^x1  ^  sin(a  +  9) 

X 

Qx+Sx 

mA+nig  COS0  — cos(q;  +  0)  cos(q!  +  0) 

y 

Qy+gy 

Ia+Ib  I  a 

9 

Qe+Se 

symm. 

a 

Qa  +  ga  . 

where  theg^.,  j  EL{x,y,6,a} ,  are  the  external  generalized  forces,  and  gj  are  the  nonlinear 

contributions  of  centrifugal,  elastic,  and  gravity  forces.  For  the  case  of  hovering  flight,  the 
motion  point  b  is  prescribed.  The  kinematics  of  {x{t),y{t),0{t))  are  prescribed  functions  of 
time.  This  reduces  the  unknowns  of  (3.1)  to  a  single  equation  fox a{t).  This  system,  which 
governs  the  evolution  of  the  deformation  of  the  profde,  has  the  form  of  a  nonlinear  oscillator; 
that  is. 


I^a  +  ka^  —1^6  -\-myq^x  sin(Q;  +  0)  +  Q„.  (3.2) 

This  equation  does  not  contain  any  structural  damping  since  no  assumptions  have  been  made 
about  a  particular  damping  model.  A  consequence  of  this  choice  is  that  this  model  cannot  be 
excited  at  (linear)  resonance.  Combes  and  Daniel  (2003a,  2003b)  proposed  that  certain  insects 
flap  their  wings  near  linear  resonance,  and  they  applied  an  arbitrary  proportional  damping  to 
make  a  finite  element  model  response  match  their  data.  The  damping  factor  used  in  the  work  of 
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Combes  and  Daniel  (2003a)  is  reported  to  be  10  times  the  mass.  If  that  same  factor  was  used  to 
impose  linear  damping  in  this  model,  the  system  would  be  overdamped.  This  level  of  damping 
would  dominate  the  response  and  the  interesting  interactions  with  the  fluid  would  be  removed. 
In  addition,  structural  damping  data  reported  in  the  literature  has  been  limited;  there  is  not  yet 
enough  evidence  to  support  a  particular  material  damping  model.  In  the  current  investigation, 
the  researchers  are  primarily  interested  in  the  elastic  response  of  the  flapping  system  with  fully 
coupled  FSI,  and  significant  damping  may  reduce  the  influence  of  aerodynamic  forces  making 
the  response  structurally  dominated.  Therefore,  structural  damping  is  not  considered  at  this  stage 
but  is  left  as  possible  addition  in  future  work. 

3.2  Fluid-structure  interactions 

The  computation  of  dynamics  of  the  fluid  field  was  carried  out  by  using  two  models  in 
conjunction  with  the  structural  model  defined  in  §3.1.  The  first  one  was  based  on  the  direct 
numerical  simulation  (DNS)  of  the  Navier-Stokes  equations  for  incompressible  flows,  and  the 
second  fluid  model  used  was  based  on  the  unsteady  vortex  lattice  method  (UVLM). 

The  coupling  scheme  is  a  predictor-corrector  method  (Preidikman,  1998;  Yang  et  ah,  2008),  as 
outlined  in  Figure  3.  Although  not  monolithic,  the  main  benefit  to  this  strategy  is  that  fluid 
model  and  structural  model  can  be  modular  and  replaced  as  needed.  The  fluid  system  is  coupled 
to  the  structural  system  by  computing  the  resultant  forces  from  the  pressure  and  vorticity  on  the 
surface  of  the  body.  This  force  is  then  used  to  integrate  a  predicted  set  of  structural  states.  In 
this  predicted  configuration,  the  surface  kinematics  of  the  body  is  fed  back  as  the  immersed 
boundary  conditions  to  the  fluid  model.  This  cycle  of  communication  is  iterated  until  the  change 
between  sub-steps  is  below  a  set  tolerance.  For  all  of  the  models  used,  only  1  or  2  sub-iterations 
are  found  to  be  needed.  This  is  likely  due  to  the  small  time  steps  used  to  comply  with  the 
Courant-Friedrichs-Lewy  condition  (CFL)  number.  This  approach  provides  a  systematic  means 
to  couple  the  equations  of  motion  in  strong  form. 
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Figure  3:  Process  flow  diagram  of  the  fluid-structure  coupling  scheme. 

3.3  Parameterization 

The  kinematics  used  are  based  on  simple  harmonic  hovering  (Wang  et  al.,  2004).  These  have 
been  adapted  to  include  a  non-impulsive  start  (Vanella  et  al.,  2009). 


x{t)  ={l  —  e  '^^)^sin(a;yt) 

y{t)  =0 

0(t)  =0Q-|-^1  —  e  j  sin((Uy. t  T  0) 


(3.3) 


Here  {x{t),y{t),6{t))  are  the  location  and  orientation  of  the  top  segment  of  the  profile  shown  in 
Figure  2.  The  exponential-decay  factor  is  used  to  prevent  the  numerical  noise  from  an  impulse 
start.  This  was  found  to  be  beneficial  in  removing  start  up  noise  in  the  simulations  and  having 
the  flow  fields  remain  well  behaved  through  the  initial  transients.  The  value  of  r  is  chosen  so 
that  it  takes  around  5  periods  of  hovering  to  achieve  the  regular  frill  motion. 
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The  nondimensionalization  of  the  system's  parameters  reduces  the  parametric  space  to  4  key 
ratios:  .  The  reference  speed  of  the  Reynolds  number  {Re  =  u^^^l ! v)  is 

Auid  ^ 

chosen  to  be  the  peak  translation  speed  of  the  forcing 


max  X  — 


The  reference  length  is  the  chord  / .  The  ratio  between  the  stroke  length  and  the  chord  length  is 
obtained  from  (Wang  et  ah,  2004)  to  be  A^l  1  =  2.8 .  The  maximum  value  for  rotation  is  set  to 

Ag=n !  A ,  and  the  profile  is  assumed  to  rotate  about  the  vertical  position  6^=-n  H .  For 
symmetric  hover  (zi  =  0 ,  and  there  is  no  lead  or  lag.  Fixing  =  1 ,  and  /  =  1  provides  a  period  of 
T  =  l.%n. 


The  density  of  the  fluid,  in  the  nondimensional  DNS  code  is  already  1;  so  to  set  the  Reynolds 
number,  the  kinematic  viscosity  v  is  selected.  The  parameters  of  the  structure  ^,k) 

Phod 

from  (3.2)  are  computed  by  geometry  and  choice  of  ( — .  In  Figure  2,  each  rigid  link  is 

Pf\mA 

taken  to  be  a  rectangle  with  a  circular  endcap.  Due  to  computational  considerations  in  the  fluid 
solver,  a  finite  thickness  profile  is  needed.  The  area,  location  of  the  centroid  ,  mass  ,  and 

rotary  moment  of  inertia  can  be  directly  computed  from  geometry  once  the  density  ratio  is 
chosen.  Finally,  the  spring  constant  k  is  computed  from  the  relation 


In  the  next  sub-sections,  the  results  obtained  through  the  DNS  and  UVLM  computations  are 
introduced  and  discussed.  For  these  computations,  the  ratio  =  25  has  been  chosen  to 

scale  the  aerodynamic  forces  to  be  of  the  same  order  as  the  fluid  forces.  The  fluid-structure 
interactions  are  investigated  by  selecting  various  spring  constants  in  the  model.  As  discussed 
above,  the  difference  spring  values  correspond  to  different  choices  of  the  ratio  o).  I  In  the 

following  results,  the  values  of  co^  /  range  from  the  soft  case  corresponding  to  the  ratio  of  1/2, 

to  the  intermediate  spring  constants  of  1/3  and  1/4,  and  the  almost  rigid  spring  case 
corresponding  to  <»^  /  <»„  =  1  /  6 . 
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3.4  Direct  numerical  simulations  in  two  dimensions 

Direction  numerical  simulation  (DNS)  represents  the  highest  fidelity  computational  fluid  model 
in  common  use.  It  is  constructed  by  the  direct  discretization  on  a  staggered  grid  of  the  Navier- 
Stokes  equations  for  incompressible  flow.  The  results,  as  presented  in  Vanella  et  al.  (2009),  are 
constructed  from  a  second  order  central  difference  scheme  on  a  stretched  Cartesian  mesh.  Time 
marching  is  performed  by  using  the  fractional  step  method  (Kim  and  Moin,  1985).  The  body  is 
represented  in  the  fixed  grid  by  using  an  immersed  boundary  method  (Yang  et  al.,  2008).  To 
enforce  the  no-slip  condition,  the  predicted  velocity  field  of  the  fractional  step  method  is  forced 
to  match  the  velocity  field  along  the  surface  of  the  body.  The  flapping  profile  is  placed  in  the 
center  of  a  large  box  so  that  the  boundaries  do  not  interact  with  the  body.  The  equations  are 
integrated  from  rest  to  15  periods  of  motion.  It  takes  approximately  1.5  days  to  compute  a  single 
period  T  of  flapping  motion  on  an  Intel  XEON  based  computer. 

The  computational  grid  was  constructed  to  resolve  the  boundary  layers,  and  other  flow  features 
on  the  moving  profile.  A  schematic  of  the  domain  is  shown  in  Figure  4  with  an  expanded  view 
near  the  tip  of  the  body.  The  center  point  of  the  profile  is  located  at  the  center  of  a  30/x30/ 
domain  to  minimize  the  effects  of  the  far-field  boundary  conditions.  The  center  region,  where 
the  body  passes  through,  is  a  uniform  grid  with  cell  size  Ax  =  Ay  =  3.725x10”^/ .  This  provides 
approximately  8  or  16  points  inside  the  boundary  for  the  various  Reynolds  number  cases. 
Outside  of  this  region,  the  grid  is  stretched  to  the  boundaries. 


Figure  4:  Schematic  of  the  profile  inside  the  DNS  grid  showing  the  overall  size  of  the  domain. 
The  detailed  view  shows  a  close  up  near  the  leading  edge  of  a  10%  thick  body,  where  the  red  lines 
are  grid  lines,  the  black  points  are  control  points  representing  the  body,  and  the  flow  field  is  the 

magnitude  of  velocity  for  Re  =  250  ,  (O^-  !  (O^  =  \  !  2  ,dX  time  t !  T  =  9.75  . 
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Figure  5:  Vorticity  contours  from  DNS  at  Re-15.  Contours  range  from  -10  to  10  with  80 
intervals.  Columns  A-C  are  flexible  profiles  with  coj' ! 1/3,  1/4}  ,  respectively;  Column 

D  is  the  rigid  profile.  Adapted  from  Vanella  et  al.  (2009). 
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Frequency  Ratio  ujf/u^n 

Figure  6:  Averaged  lift  to  drag  ratio  from  the  two-dimensional  DNS  calculations  for  various 
frequency  ratios.  Profile  thickness  is  I/IO. 

A  summary  of  the  resulting  flow  fields  is  depicted  in  Figure  5  for  a  range  of  /  0)„  at  Reynolds 

number  Re  =  75 .  The  vorticity  contours  reveal  the  vortex  structure  interplay  with  the  flexible 
profile.  Similar  results  were  also  computed  for  Re  =  250  and  Re  =  1000  .  It  is  worth  noting  that 
for  the  soft  spring  case  of  =  1  /  2 ,  the  system  undergoes  extremely  large  deflections,  and 

the  passive  link  almost  undergoes  a  complete  rotation  about  the  joint. 

On  examining  the  averaged  dimensionless  aerodynamic  forces  acting  on  the  profile,  there  is  an 
interesting  finding.  It  is  observed  that  for  the  particular  spring  value  corresponding  to 
&)y  /  &)„  =  1  /  3 ,  there  is  a  peak  in  the  ratio  of  lift  coefficient  and  drag  coefficient  Q  /  ,  as 

viewed  in  Figure  6.  It  is  also  noted  that  the  flexible  profile  has  an  improved  efficiency  compared 
to  that  of  the  rigid  profile. 

3.5  Unsteady  vortex  lattice  method  in  two  dimensions 

In  contrast  to  the  computationally  expensive  DNS  method,  vortex  methods  present  a  compromise 
between  speed  and  fidelity.  In  the  unsteady  vortex  lattice  method  (UVLM)  employed  by 
Preidikman  (1998),  it  is  assumed  that  the  flow  field  is  inviscid  and  the  wake  can  be  completely 
described  by  point  vortices.  A  body  in  the  flow  field  is  discretized  into  panels,  and  the  no¬ 
penetration  condition  is  applied  at  the  chosen  control  points  along  each  panel.  The 
discretization  used  in  this  method  is  illustrated  in  Figure  7. 
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Figure  7:  Schematic  for  the  discretization  of  the  Unsteady  Vortex  Lattice  Method. 

This  method  can  be  applied  to  a  membrane  or  zero-thickness  problems,  since  assuming  the  wake 
to  only  separate  at  the  edges  is  suitable.  Vortices  are  convected  from  the  trailing  edge  after  each 
time  step  by  using  the  Kutta  condition.  Similarly,  Valdez  et  al.  (2006)  proposed  a  method  to 
convect  vortices  simultaneously  from  the  leading  of  the  profile.  Since  a  zero-thickness  profile  is 
used  here,  the  tips  are  assumed  to  be  the  points  where  the  wake  separates  from  the  body.  Also 
proposed  by  Valdez  et  al.  (2006),  is  a  reconstruction  of  the  entire  velocity  and  pressure  fields 
from  the  vortex  particle  wake.  For  the  results  presented  here,  the  original  code  was  authored  by 
Valdez  (2008). 

Unlike  the  slower  DNS  computations,  a  complete  run  of  20  periods  of  hovering  motion  takes 
around  10  hours  on  an  eight  processor  Intel  XEON  computer.  The  issue  with  longer  simulations 
is  that  with  each  additional  time  step,  there  are  two  additional  vortices  whose  influence  needs  to 
be  included.  So  the  number  of  computations  increases  at  a  rate  0{n^) ,  where  n  is  the  number  of 
time  steps  in  the  integration.  This  makes  short  time  computations  quick  in  comparison  to  a  DNS 
study,  but  long  time  simulations  quickly  become  impractical.  These  characteristics  make  the 
UVLM  simulations  attractive  from  a  design  perspective  since  coarse  results  can  be  obtained 
rather  quickly. 

3.6  POD  analysis  of  the  flow  fields 
3.6.1  Formulation 

The  proper  orthogonal  decomposition  (POD)  goes  by  many  names  such  as  the  Karhunen-Loeve 
transform,  principal  component  analysis,  or  singular  systems  analysis  depending  on  the 
discipline.  It  also  can  be  formulated  in  a  continuous  or  discrete  sense,  and  used  for  experimental 
or  computational  data.  The  overarching  goal  of  the  POD  is  to  decompose  data  into  hierarchical 
sets  of  spatial  basis  functions,  often  called  mode  shapes.  Here  the  velocity  fields  of  the  fluid  are 
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decomposed  into  spatial  modes  by  using  the  continuous  approach.  The  gives  the  assumed 
representation  of  the  flow  field  the  following  form 

oo 

=  0,(x).  (3.4) 

i=l 

This  can  be  interpreted  from  a  vibrations  perspective  as  time-dependent  coefficients  in  modal 
coordinates.  The  modes  ^(x)  are  chosen  to  maximize  the  projection  of  the  empirical  data  onto 

these  modes  in  a  sense.  So  the  problem  statement  is  formulated  as  in  the  work  of  Holmes  et 
al.  (1996),  and  further  simplifications  can  be  made,  since  it  can  be  recognized  that  the  modes  are 
a  special  superposition  of  the  data  snapshots.  Employing  the  method  of  snapshots  (Sirovich, 
1987),  the  modes  (j)  can  be  approximated  as  a  finite  sum  over  the  known  data 


M 


k=\ 


m(x,4) 

v(x,4) 


(3.5) 


O’Donnell  and  Helenbrook  (2007)  demonstrate  through  a  scaling  argument  that  the  pressure 
components  can  be  neglected  for  incompressible  flow  in  the  substitution  back  into  the  Fredholm 
equation.  This  means  that  only  the  velocity  field  is  needed  for  the  computations,  while  the 
modes  of  the  pressure  field  are  still  computed.  The  simplified  results  become  an  algebraic 
eigenvalue  problem  of  size  M 


Cf  =  Xf 


(3.6) 


where 


Qi  :=^  f  «(x,6)- u(x,4)fi?r.  (3.7) 

M 

Since  C  =  only  the  upper  or  lower  triangular  part  needs  to  be  constructed.  So  the  process  to 
construct  the  POD  set  can  be  described  as  a  sequence  of  the  following  steps: 

•  Generate  velocity  field  data  at  equal  time  intervals. 

•  Construct  each  element  of  (3.7)  by  integrating  over  the  domain. 

•  Solve  the  algebraic  eigenvalue  problem  of  (3.6)to  get  the  set  of  eigenvalues  X  and 

the  associated  eigenvectors  ip . 

•  Back  substitute  y/  into  (3.5)  to  get  a  truncated  set  of  POD  eigenfunctions  (j) . 
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The  set  of  spatial  mode  shapes  (j)  are  ranked  by  using  the  eigenvalues  X .  This  gives  a 
quantitative  measure  to  regard  the  importance  of  each  mode  as  constructed  from  the  data.  For  all 
the  examples  considered  below,  99%  of  the  energy  in  the  snap  shots  in  contained  in  less  than  10 
of  the  leading  eigenvalues.  This  significant  roll-off  means  that  one  can  consider  just  a  few 
leading  terms.  A  noteworthy  side  effect  of  the  POD  obtained  by  using  the  method  of  snapshots 
is  immediately  appreciable;  (j)  must  be  divergence  free.  Since  the  mode  shapes  are  a  weighted 
super  position  of  divergence  free  data,  then  (j)  must  also  be  divergence  free  (incompressible).  So 
any  use  of  these  modes  faithfully  preserves  the  incompressibility  of  the  flow  field. 

3.6.2  Results 

Next,  the  results  obtained  from  various  POD  computations  are  presented.  Note  that  these  figures 
represent  the  flow  field  outside  of  time,  and  they  are  the  hierarchical  structures  of  the  fluid  flow. 
The  center  joint  of  the  profile  moves  through  the  region  (x  /  l,y  1 1)  =  (+1.4,0) . 

In  Figure  8,  the  vorticity  contours  of  the  computed  modes  are  shown  fori?e  =  75.  The  plotted 
field  is  the  curl  of  the  velocity  mode  shape,  normalized  to  have  the  maximal  value  be  one.  This 
makes  comparison  between  modes  informative.  After  examining  this  figure,  it  is  observed  that 
all  of  the  mode  one  results  contain  a  large  pair  of  vortices.  This  represents  the  downward  jet  of 
fluid,  and  produces  lift  on  the  profile.  For  =  1  /  3 ,  this  pair  of  vortices  is  the  most  intense 

as  well  as  closest  to  the  region  of  the  body.  Interestingly  this  corresponds  to  the  most  efficient 
response  for  the  given  harmonic  input  kinematics,  ft  also  corresponds  to  the  structural  resonance 
for  the  kinematics  being  used.  As  expected,  the  scale  of  the  structures  decreases  as  the  mode 
number  increases.  This  follows  intuition  from  vibration  mode  shapes. 


13 

Approved  for  public  release;  distribution  unlimited. 


-2  -1  0  1  2-2-10  1  2-2  -1  (1  1  2 

-1  -0.8  -n.fi  -0.4  -0.2  1)  0.2  0.4  O.fi  0.8  i 

Figure  8:  Vorticity  contours  determined  by  POD  from  DNS  at  Re  =  15 , 10%  thickness,  and 
harmonic  kinematics.  Normalization  is  maxi  V  x  (p:  (x)  1=  1 .  Florizontal  scale  is  x !  I ,  vertical 

X  ' 

scale  is  y !  I . 

The  results  obtained  for  Re  =  250  are  shown  in  Figure  9.  It  is  interesting  to  note  similar  patterns 
to  those  obtained  for  the  low  Reynolds  number  eonfiguration;  again,  mode  one  eomponents 
shows  a  downward  jet  and  mode  two  shows  the  end  of  stroke  vortex  pair.  In  this  flow  regime, 

the  vortiees  for  the  spring  of  intermediate  stiffness  ( e  {1  /  3, 1  /  4} )  appears  to  be  more 

spatially  regular.  An  investigation  using  other  kinematics  would  need  to  be  carried  out  in  order 
to  see  how  the  kinematics  affects  the  regularity  of  the  field. 
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Figure  9:  Vorticity  contours  determined  by  POD  from  DNS  at  Re  =  250 , 10%  thickness,  and 
harmonic  kinematics.  Normalization  is  maxi  V  x  (p:  (x)  1=  1 .  Florizontal  scale  is  x !  I ,  vertical 

X  ' 

scale  is  y !  I . 

Comparing  the  mode  shapes  of  DNS  fields  to  the  UVLM  fields  cannot  be  done  with  vorticity 
contours  since  the  UVLM  approach  is  inviscid  in  nature.  Instead  the  velocity  contours  must  be 
compared.  As  discussed  later  in  §3.7,  the  UVLM  velocity  modes  share  similar  features  to  the 
DNS  generated  modes. 
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3.7  Comparing  the  results  from  DNS  and  UVLM  studies 

The  discussion  in  this  section  follows  from  the  published  work  of  Fitzgerald  et  al.  (2011).  The 
goal  is  to  compare  the  DNS  and  UVLM  models  to  see  where  and  how  each  should  be  best 
employed.  A  DNS  simulation  is  expensive;  in  fact,  it  would  take  several  weeks  for  a  full  run  of 
15  periods.  By  contrast,  the  UVLM  can  produce  15  periods  of  flapping  in  several  hours.  The 
DNS  results  are  expected  to  be  better,  since  it  fully  models  the  physics  of  interest.  However,  as 
shown  by  Fitzgerald  et  al.  (2011)  the  trends  observed  in  the  UVLM  simulations  for  various 
configurations  agree  well  with  those  noted  in  the  DNS  studies.  This  provides  a  hybrid  approach 
to  use  for  system  designers:  the  UVLM  simulations  can  be  first  used  as  a  low-fidelity  prediction 
tool  to  find  parametric  regions  of  interest,  which  can  be  followed  by  the  use  of  DNS  studies  to 
compute  more  realistic  data. 

3.7.1  Flow  fields 

In  Figure  10,  comparisons  among  the  flow  fields  obtained  through  the  DNS  studies  for  Re  =  15, 
250  and  1,000  and  the  UVLM  are  shown  for  period  6.  The  ratio  of  the  forcing  frequency  to  the 
natural  frequency  of  the  system  is  chosen  to  be  <x>j  !  co^=\l  2  to  show  the  largest  displacements 

observed.  The  magnitude  of  the  velocity  field  is  shown,  since  the  vorticity  field  is  unavailable 
from  the  UVLM  simulations.  In  these  figures,  the  different  vortex  structures  can  be  observed  and 
compared  as  the  viscous  diffusion  in  the  system  is  decreased. 
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Figure  10:  Magnitude  of  the  velocity  fields  from  periods  6—7  for  co^  !  =  \ll ,  DNS  at 

Re  e  {75, 250, 1000} ,  and  UVLM.  The  velocity  field  is  normalized  as  |  u  . 

For  Re  =  15,  the  flow  field  snapshots  throughout  periods  6  and  9  are  nearly  identical,  resulting  in 
the  periodicity  of  the  flow  field.  This  matches  with  the  calculations  of  the  correlation  dimension, 
which  suggests  that  a  periodic  orbit  is  produced  in  the  low  Reynolds  number  case.  When 
Re  =  250,  the  snapshots  are  also  very  similar  between  periods  6  and  9.  However,  as  the 
correlation  dimension  indicates,  an  exact  periodic  solution  is  not  seen.  There  are  small  enough 
differences  in  the  flow  fields  to  cause  small  disturbances  to  what  is  nearly  a  periodic  orbit.  At 
Re  =  1000 ,  with  relatively  low  viscous  diffusion,  the  flow  field  no  longer  appears  to  be  periodic. 
However,  the  same  kind  of  vortex  structures  can  be  identified  as  in  the  low  Reynolds  number 
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cases,  such  as  the  leading  and  trailing  edge  vortices.  As  expected,  the  intensity  of  these  vortex 
structures  increases,  since  they  are  proportional  to  the  magnitude  of  the  velocity.  The  non¬ 
periodicity  of  the  flow  field  is  a  consequence  of  vortex  interactions  that  were  not  dominant  in  the 
low  Reynolds  numbers  cases. 

A  sample  of  the  velocity  POD  mode  shapes  obtained  from  the  DNS  and  UVLM  studies  is  shown 
in  Figure  11.  Here,  the  downward  jet  appears  in  the  vertical  velocity.  Around  the  body,  traces  of 
the  end  of  stroke  vortices  can  be  noted.  The  UVLM  contains  no  viscous  dissipation  which  is 
why  the  fluid  structures  do  not  dissipate  in  time.  The  results  show  that  the  UVLM  POD  is  still 
able  to  nicely  quantify  the  flow  field  in  terms  of  the  dominant  structures. 


He  =  75 


He  =  250 


UVLM 
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Figure  11:  Comparisons  of  Mode  1  of  POD  velocity  contours  from  DNS  and  UVLM  data  for 
co^  !  co^  =1/3  with  harmonic  kinematics.  The  fields  are  normalized  by  max  |  (j).  (x)  1^  =  1 . 

X 

Horizontal  scale  is  x  /  / ,  vertical  scale  is  y !  I . 
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3.7.2  Aerodynamic  loads 

The  aerodynamic  forces  of  interest  here  are  the  lift  force  and  the  drag.  The  lift  is  defined  as  the 
vertical  force;  or  the  force  that  would  keep  a  hovering  system  from  falling.  The  dimensionless 
lift  is 


Q 

Q  =  ^  (3.8) 

Pfluid^ref  ^ 

where  is  the  vertical  force.  Since  the  vertical  motion  of  the  joint  is  zero,  the  drag  is  taken  to 
be  the  horizontal  force  that  opposes  the  horizontal  translation.  The  dimensionless  drag  is 

C^=-sgn(i)^  (3.9) 

^  Pfluid^ref  ^ 

where  is  the  horizontal  force.  Observing  the  time  series  obtained  through  DNS  at  Re  =  1000 

and  the  UVLM  shows  that  the  signals  initially  match  quite  well,  as  seen  in  the  first  period  and  a 
half  in  Figure  12.  As  the  simulation  evolves  however,  the  solutions  appear  to  diverge,  but  still 
follow  similar  trends.  The  loads  from  the  UVLM  successfully  capture  the  dominant 
characteristics,  with  similar  peaks  and  valleys  in  Q(0  corresponding  to  the  mid-stroke  and 

stroke  reversals,  respectively.  In  both  the  lift  and  drag  coefficients  the  UVLM  results  appears  to 
follow  the  same  pattern  of  over-estimating  the  highs  and  lows.  This  suggests  that  the  simulation 
of  a  transient  maneuver,  such  as  clap  and  fling  (Weis-Fogh,  1973),  and  averaged  long-term 
dynamics  may  be  reasonably  predicted  by  using  UVLM  simulations.  To  explore  long-term 
dynamics,  a  statistical  approach  is  needed. 
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Figure  12:  Time  series  of  the  dimensionless  loads  of  the  hovering  profile  from  DNS  at  Re  =  1000 
in  red  and  UVLM  in  blue  for  the  frequency  ratio  co^  !  co^  =  1  /  3 . 

The  time  histories  have  been  collapsed  to  a  single  period  by  phase-averaging  in  Figure  13.  The 
data  from  periods  5  through  15  are  averaged  by  using  the  hovering  period  T  as  the  reference 
clock.  This  range  of  time  is  used  since  the  hovering  kinematics  has  reached  the  full  amplitude, 
and  the  startup  transients  in  the  fluid  should  have  died  down.  At  the  stroke  reversal,  when 
t !  T  e  (0,  1  /  2} ,  there  is  a  jump  in  that  is  not  fiilly  captured  by  the  UVLM  simulations.  A 

possible  reason  for  the  discrepancy  could  be  that  the  vortex  interaction  during  this  event  is 
influenced  more  by  viscosity  than  at  other  points  of  the  cycle.  For  softer  spring  values;  that  is, 
e  {1  /  2,  1  /  3,  1  /  4} ,  the  curves  generally  are  in  better  agreement.  This  is  particularly 

visible  in  the  case  of  the  lift.  The  prediction  is  worst  for  the  rigid  case. 
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Figure  13:  Phase  averaged  forces  of  the  hovering  profile  from  DNS  at  Re  =  1000  in  red  and 
UVLM  in  blue  for  various  frequency  ratios,  (a)  Lift  coefficient  and  (b)  drag  coefficient.  Adapted 

from  Fitzgerald  et  al.  (2011). 

The  quantitative  comparison  of  which  UVLM  curve  best  matches  the  corresponding  DNS  curve 
requires  another  definition.  The  error  measure  assumed  here  is  a  scaled  -norm.  This  is  a 
point-wise  check  to  see  how  far  off  the  UVLM  is  at  every  point  in  time.  This  does  not  account 
for  phase  lag  or  lead  between  the  models.  For  the  drag  force,  it  is  written  as 


error  (Cp) 


(3.10) 


with  the  same  form  being  used  for  the  lift  coefficient  Q .  The  numerical  values  are  compiled  in 


Table  1.  The  differences  appear  to  be  smaller  for  the  more  compliant  structures,  indicating  that 
the  UVLM  model  would  be  of  better  use  in  highly  flexible  configurations  where  the  structural 
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Table  1 :  Relative  error  values  of  the  phase-averaged  loads  between  the  Re  =  1 000  DNS  and 

UVLM  results  for  a  range  of  stiffiiesses. 


Frequency  Ratio  ujj-  / 


Rel.  Error 

1/2 

1/3 

1/4 

1/6 

Rigid 

c,, 

9.9% 

15.9% 

10.0% 

19.6% 

24.7% 

C, 

9.6% 

13.4% 

13.6% 

25.6% 

28.8% 

-•-Clint  =  7^  ■  Ci,;  7?e  =  250  Rt  =  1000  UVLM 


Figure  14:  Comparisons  of  time  averaged  dimensionless  lift  and  drag  coefficients  from  the  UVLM 
and  DNS  at  various  Re .  (a)  Mean  lift  and  drag  coefficients,  (b)  the  ratio  of  mean  lift  to  mean 
drag.  Adapted  from  Fitzgerald  et  al.  (201 1). 

dynamics  outweigh  the  fluid  contributions.  Overall,  the  predictions  are  likely  acceptable  in 
engineering  design  since  the  compliant  configuration  have  errors  that  are  less  than  20%. 

Taking  the  overall  time-average  of  the  loads  provides  an  even  better  use  of  the  UVLM  results. 
As  shown  in  Figure  13,  the  time  averaged  loads  Q  and  trend  together.  The  UVLM 

simulations  overestimate  the  lift  and  the  drag,  but  since  it  is  nearly  by  the  same  15%  the  ratio  of 
the  quantities  match  the  DNS  results  quite  well.  Keeping  the  modeling  and  kinematic 
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assumptions  in  mind,  the  results  indicate  that  moderate  flexibility  improves  the  aerodynamic 
performance. 

The  results  of  Figure  13  indicate  that  the  UVLM  study  is  adequate  for  quasi-steady  prediction. 
Furthermore,  the  results  of  Figure  14b  demonstrate  that  the  UVLM  is  quite  suitable  for  a 
designer  to  predict  gross  quantities  and  trends.  The  trends  with  respect  to  the  spring  parameter 
are  the  most  significant,  since  it  points  to  the  UVLM  study  being  suitable  for  use  in  an 
optimization  setting.  As  mentioned  earlier,  a  designer  could  adopt  a  hybrid  approach  and  use  the 
UVLM  simulations  in  an  optimizer  to  find  parametric  regions  of  interest  and  follow  this  up  with 
detailed  investigations  of  these  regions  with  DNS  studies. 


23 

Approved  for  public  release;  distribution  unlimited. 


4  Three-Dimensional  Modeling 

4.1  Introduction 

A  survey  of  the  technology  to  extend  the  flapping  work  into  3D  problems  quickly  reveals  that  the 
modeling  of  the  structure  is  less  well  understood  than  the  modeling  of  the  fluid.  The  Navier- 
Stokes  equation  for  an  incompressible  Newtonian  fluid  is  an  extremely  versatile  and  well- 
understood  model.  Even  a  vortex  lattice  method  in  3D  is  well  suited  for  certain  problems.  The 
model  for  a  body  or  structure,  by  contrast,  is  still  an  open  item.  Not  merely  the  geometry  is 
variable,  but  the  structure  itself  and  how  to  model  its  behavior  is  open.  Recent  work  on 
modeling  the  flapping  wing-fluid  system  have  either  included  rigid  wings  or  highly  simplified 
models.  In  the  work  of  Pai  et  al.  (2009),  a  body  described  by  nonlinear  structural  elements 
interacts  with  a  quasi-steady  fluid.  In  an  experimental  study,  Zhao  et  al.  (2010)  measured  the 
forces  on  a  flapping  wing  made  with  Mylar  material.  The  synthesis  of  natural  fliers  to  the 
construction  of  micro  aerial  vehicles  is  also  a  topic  of  open  interest. 

The  goal  of  this  section  is  to  present  a  structural  model  that  is  modular  enough  to  be  adapted  to  a 
variety  of  geometry  and  material  properties,  and  can  be  designed  for  use  in  concert  with  a  large- 
scale  fluid  simulation.  The  finite  element  method  (FEM)  provides  a  nice  foundation  to  base  the 
model  since  there  is  extensive  literature  on  many  aspects  of  its  use.  The  FEM  is  also  adaptable 
to  nearly  any  geometry  and  has  the  potential  to  be  generalized  for  a  wide  range  of  material 
systems. 

For  the  immersed  boundary  methods  employed  in  the  available  large-scale  fluid  solver,  a  body 
with  finite  thickness  is  required.  It  was  this  same  requirement  that  necessitated  the  placement  of 
a  surface  around  the  structural  elements  in  the  2D  work  of  the  previous  section.  Several  fluid 
grid  points  must  be  considered  inside  the  solid  body  for  the  pressure  to  be  resolved  properly. 
Therefore,  the  body  description  selected  was  a  solid  body,  and  not  a  plate  or  shell.  This  provides 
a  natural  thickness  to  the  model,  and  allows  for  future  cases  with  highly  detailed  surface 
geometry  such  as  a  computed  tomography  scan  of  an  insect  wing.  A  solid  element  also  has  the 
benefit  that  one  can  directly  employ  many  different  material  laws  from  continuum  mechanics. 

4.2  Experiments 

The  bio-inspired  design  of  flapping  vehicles  draws  heavily  from  the  work  of  both  engineers  and 
biologists.  Among  the  various  interesting  aspects  of  insect,  the  wing  is  one  of  them.  Largely 
thought  to  be  a  passively  flexible  structure,  the  details  of  the  structure  of  the  wing  have  long  been 
an  area  of  interest.  Comstock  (1918)  has  illustrated  a  life's  body  of  work,  which  includes 
cataloging  and  defining  of  the  various  characteristics  of  insect  wing  morphology.  His  naming 
conventions  of  the  venation  are  still  in  common  use  by  biologists  today.  Following  these 
footsteps,  the  efforts  of  Wootton  (1992)  further  built  on  the  biological  map  of  insect  wings. 
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Early  engineering-type  studies  to  model  the  wing  include  the  efforts  of  Combes  and  Daniel 
(2003a,  2003b,  2003c)  and  the  PhD  dissertation  of  Wootton's  student  (Herbert,  2001).  Both 
groups  come  from  a  biology  background,  and  they  make  claims  that  the  flapping  frequency 
observed  is  the  fundamental  frequency  of  the  structure  of  the  wing.  This  notation  that  an  insect 
flaps  at  linear  resonance  has  only  recently  been  challenged.  Computationally,  this  does  not  agree 
with  the  results  shown  in  §3.7.2  since  those  models  predict  that  efficiency  decreases  as  linear 
resonance  is  approached  (Fitzgerald  et  ah,  2011;  Vanella  et  al.,  2009). 

Experimentally,  Sims  (2010)  showed  that  the  first  natural  frequency  of  a  hawkmoth  wing  is 
around  twice  that  of  the  flapping  frequency.  The  spectral  information  was  measured  via  a 
scanning  laser  vibrometer  from  wings  recently  removed  from  living  hawkmoths.  The  tests  were 
repeated  for  several  specimens  in  air  and  in  a  vacuum  chamber,  and  the  measured  first  natural 
frequency  is  around  twice  the  flapping  frequency.  Kang  et  al.  (2011)  made  several  scaling 
arguments  that  predict  the  optimal  natural  frequency  to  be  around  1/3  to  1/2  of  the  flapping 
frequency.  This  recent  work  indicates  that  not  only  could  older  assumption  be  wrong,  but  there 
is  a  tremendous  opportunity  to  exploit  nonlinear  effects. 

The  setup  to  measure  the  spectral  response  of  a  living  insect  is  outlined  in  Figure  15.  A  living 
insect  is  anesthetized  by  exposure  to  a  large  amount  of  Flynap^  and  then  placed  in  a  fixture 
molded  out  of  modeler's  clay.  The  forewing  is  fixed  in  place  near  the  root  with  more  clay  and 
pins.  The  wing  is  acoustically  excited  by  a  JBL  ASB1728  loudspeaker.  This  subwoofer  is  rated 
to  4000  W  of  continuous  pink  noise,  and  has  high  fidelity  down  to  20  Hz.  A  pseudo-random 
signal  is  output  from  the  Polytec  PSV-400  vibrometer  controller  to  a  Crown  MA-9000i  power 
amplifier  connected  to  the  loudspeaker. 

The  choice  of  the  Manduca  sexta  was  made  due  to  several  key  factors.  It  has  been  widely 
studied  by  biologist  and  found  to  be  rather  uniform  in  its  body  and  flight  characteristics  across 
individuals.  The  insect  wing  is  relatively  large  and  opaque.  This  means  that  it  can  measured  by 
standard  vibrometer  equipment  already  available  in  the  Vibrations  Laboratory.  This  insect 
species  can  be  easily  procured  and  grown  from  larvae  purchased  from  biology  supply  companies. 
During  the  design  of  the  experiment,  it  was  thought  that  the  fundamental  frequency  of  these 
insects  was  near  25  Hz  (Combes  and  Daniel,  2003a,  2003b,  2003c).  This  drove  the  interest  in 
the  use  of  a  speaker  with  good  low  frequency  fidelity. 


^  Flynap  is  a  general  anesthetic  designed  for  Musca  domestica.  It  is  composed  of  50%  Triethylamine,  25%  Ethanol, 
and  25%  Fragrance,  per  the  Material  Safety  Data  Sheet  httr)://www.carolina.com/r)df/msds/FLYNAP.pdf 
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Loudspeaker 

Figure  15:  Schematic  depicting  the  use  of  a  scanning  laser  vibrometer  to  characterize  the  spectral 
response  of  a  living  insect  wing  using  non-contact  excitation. 


Figure  16:  Mesh  of  the  scanning  laser  vibrometer  on  a  living  Manduca  sexta  forewing.  The  color 
markers  indicate  locations  of  signal  points  in  Figure  17.  Note  that  +x  is  in  the  vertical  direction. 

The  scanning  laser  vibrometer  is  then  setup  to  measure  the  response  of  a  set  of  points  on  the 
surface  of  the  wing.  This  mesh  is  shown  in  Figure  16.  The  laser  vibrometer  is  placed  such  that 
the  wing  is  centered  and  parallel  in  the  viewfinder.  The  scanning  is  performed  sequentially,  and 
the  frequency  information  from  each  point  is  stored  as  complex  Fast  Fourier  Transform  (FFT) 
data.  Since  the  total  number  of  FFT  data  points  per  mesh  point  is  limited  by  the  software,  there 
is  a  trade-off  of  spectral  resolution  when  selecting  the  frequency  range  of  interest.  By  choosing 
the  maximum  number  of  FFT  data  points  at  6400,  and  selecting  the  frequency  range  of  interest  to 
be  0-1000  Hz,  a  working  resolution  around  0.25  Hz  was  obtained. 

By  examining  the  FFT  data  from  several  key  points  around  the  wing,  the  first  natural  frequency 
of  the  wing  can  be  easily  located.  In  Figure  17,  the  normalized  FFT  results  are  shown  for  points 
near  the  root,  the  tip,  and  the  trailing  edge.  The  measurements  made  at  these  spatial  points 
indicate  that  the  first  natural  frequency  of  the  wing  specimen  is  around  77  Hz.  Locating  the 
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Figure  17:  Normalized  magnitude  of  the  FFT  of  the  velocity  data  determined  from  the  laser 
vibrometer  at  various  points  on  the  forewing  of  a  Manduca  sexta.  The  color  points  are  indicated 

on  Figure  16. 

second  natural  frequency  is  a  bit  more  challenging,  since  the  response  of  the  wing  appears  is 
found  to  be  rather  complicated. 

By  searching  through  the  visualizations  of  the  stitched  mode  shapes  on  the  PSV-400  system,  it  is 
found  that  around  134  Hz,  there  is  another  mode  shape.  The  noise  floor  is  too  high  at  the  higher 
frequencies  to  be  confident  in  locating  other  natural  frequencies  and  associated  mode  shapes. 
This  problem  is  inherent  in  the  non-contact  excitation  provided  by  using  a  speaker.  The  natural 
frequencies  reported  here  agree  well  with  those  of  Sims  (2010),  who  used  an  amputated  wing 
directly  mounted  to  a  shaker  and  could  therefore  work  with  much  higher  frequencies  and 
amplitudes  of  excitation.  The  novelty  of  the  speaker  experiments  is  that  the  insect  and  wing 
were  alive  during  and  after  the  entire  measurement  process.  The  wing  measured  here  has  not 
been  altered  by  death,  atrophy,  or  temperature.  Sims  had  to  measure  the  severed  wings  within 
several  hours  to  ensure  that  they  had  not  changed  significantly,  whereas  in  the  current  study,  the 
researchers  had  around  30  minutes,  the  window  in  which  the  insect  remained  asleep. 

Reconstructions  of  the  mode  shapes  associated  with  the  first  two  natural  frequencies  are  shown 
in  Figure  18.  Here,  the  data  was  extracted  from  the  proprietary  Polytec  data  file,  and  plotted  in 
Matlab.  The x and y coordinates  are  scaled  by  the  span  of  the  wing/,  and  the  vertical 
displacement  of  the  mesh  is  scaled  such  that  the  maximum  is  /  /  8  .  The  choice  of  vertical  scaling 
is  arbitrary  since  the  representation  is  for  a  mode  shape,  and  the  choice  of  /  /  8  was  merely  made 
for  visualization  purposes.  The  first  mode,  which  is  shown  in  Figure  18a,  appears  to  correspond 
to  spanwise  bending.  The  second  mode,  which  is  shown  in  Figure  18b,  appears  to  be  a 
combination  of  chordwise  bending  and  some  bending  near  the  tip.  A  possible  use  for  this  type  of 
detailed  information  is  in  the  construction  of  wing  models  tailored  to  perform  like  a  Manduca 
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Figure  18:  Experimentally  determined  modes  of  a  living  Manduca  sexta  forewing.  The  x  and  y 

coordinates  are  scaled  by  the  wing  span,  z  is  scaled  such  that  the  max  z !  I  =  \l%  for 
visualization.  (A)  Mode  1  at  77.5  Hz.  (B)  Mode  2  at  133.75  Hz. 

sexta.  The  distribution  of  material  properties  could  be  designed  such  that  the  first  two  natural 
frequencies  of  the  model  resemble  the  experimental  results. 

4.3  Implementation  using  geometrically  exact  finite  elements 

The  goal  of  the  tools  constructed  here  are  to  explore  fluid-structure  interaction  problems,  like 
flapping  wings.  Solid  finite  elements  are  used  since  the  supporting  theory  and  technology  is 
widely  known.  They  provide  a  foundation  to  build  a  framework  that  many  different  types  of 
structures,  wings,  material  models,  and  so  on,  can  be  tested.  The  implementation  is  general 
enough  to  handle  structural  elements  as  well,  and  their  integration  is  a  possible  avenue  for  future 
work.  Solid  models  were  chosen  since  they  provide  finite  thickness,  fundamentally  have  the 
fewest  assumptions,  and  they  can  be  widely  adapted  for  a  variety  of  continuum-mechanics  based 
response  models.  The  relative  cost  of  using  solid  elements  in  a  large-scale  CFD  simulation  is  not 
low. 

The  implementation  technology  employed  here  is  largely  based  on  Hughes  (2000)  for  the 
assembly  and  shape  functions  and  Belytschko  et  al.  (2000)  for  dealing  with  nonlinear  models. 
Mesh  generation  is  designed  around  the  open-source  software  Gmsh  (Geuzaine  and  Remade, 
2009).  The  elements  implemented  are  isoparametric  quadratic  hexahedra  for  the  volume  of  the 
body.  These  27-node  displacement  based  elements  were  selected  since  they  do  not  suffer  from 
locking  like  linear  elements.  Surface  elements  used  for  the  FSI  and  other  loading  are  9  node 
quadrilaterals.  Each  quadrilateral  is  coincident  to  a  single  face  of  a  corresponding  hexahedron. 

The  implementation  contains  both  Kirchhoff  and  Biot  material  models  for  isotropic  properties. 
The  implementation  of  the  Biot  material  was  selected  since  it  represents  a  generalization  of  the 
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engineering  stress-strain  law  to  finite  deformations  with  arbitrary  rigid  body  motions.  This 
method  was  selected  since  it  provides  a  pattern  for  future  code  development  for  problems  with 
anisotropic  materials. 

4.3.1  Equations  of  motion 

The  description  of  motion  implemented  is  a  weak-form  of  momentum  conservation  often  called 
the  Total  Lagrangian  formulation.  It  is  a  Lagrangian  method  where  everything  is  expressed  in 
terms  of  the  reference  configuration.  In  the  usual  finite  element  way,  the  virtual  work  of  the 
body  is  expressed  as 


=0  ■  (4.1) 

The  virtual  internal  work  can  be  expressed  in  terms  of  any  work-conjugate  pair.  The  simplest 
pair  to  work  is  (E,S) ,  which  in  Voigt  notation  takes  the  form 

=  f  {Sef{S}dV,=^6W:,.  (4.2) 

The  external  work  can  be  viewed  as  the  sum  of  work  due  to  body  forces,  such  as  acceleration 
and  gravity,  and  surface  tractions. 


(4.3) 

The  surface  forces  will  be  treated  in  §4.3.2.  The  acceleration  term  can  be  written  as 

SW,  =  f Jui-Up„yF„.  (4.4) 

By  using  the  usual  finite  element  shape  function  approximations 


u  =  Nq 


(4.5) 


one  obtains  the  consistent  mass  matrix. 

=  f  /!q’'NU-NqpJdF„  = 

M=  f 

Notable  features  of  this  matrix  are  that  it  is  symmetric,  and  constant.  The  symmetry  allows  for 
computational  efficiency  in  terms  of  storage  and  inversion.  The  fact  that  it  is  constant  means  it 


(4.6) 
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only  needs  to  be  built  once.  In  explicit  dynamic  algorithms,  it  needs  to  only  be  decomposed 
once  as  well.  This  formulation  is  still  consistent  with  large  motions,  and  no  implied  assumptions 
about  the  body  have  been  made.  The  principle  of  the  conservation  of  mass  can  be  used  to  show 
that  this  constant  mass  matrix  in  the  Total  Lagrangian  formulation  is  equivalent  to  the 
deformation  dependent  mass  matrix  by  other  formulations  like  the  updated  Lagrangian  form 
(Belytschko  et  ah,  2000). 

After  combining  (4.2)  and  (4.6)  into  (4.1),  the  result  is 

0  =  «qyf,„,+Mq-f„). 


After  inclusion  of  linear  damping,  this  can  be  recast  as  the  semi-discrete  equation  of  motion 

Mq(0  +  Dq(0  +  (q,  t)  =  f,,,  (q,  t)  (4.7) 

The  calculation  of  f;^(  is  determined  by  how  the  selected  material  model  is  used  to  compute  the 

second  Piola-Kirchhoff  stress  S .  Essential  boundary  conditions  have  not  yet  been  applied  to  this 
equation  and  these  boundary  conditions  are  needed,  along  with  initial  conditions,  to  fully  pose 
the  problem. 

4.3.2  Application  of  essential  boundary  conditions 

The  degrees  of  freedom  (DOF)  of  the  entire  body  are  ordered  during  the  preprocessing  of  the 
mesh  to  place  the  restrained  components  at  the  end  of  the  global  list.  Thus  if  q  are  the  DOF  for 
the  entire  body,  then,  this  list  is  partitioned  as 


q  = 


q 

V 


(4.8) 


where  the  q  are  the  unrestrained  DOF,  and  v  are  the  DOF  that  have  some  essential  boundary 
condition  applied  to  them.  Here,  \(t)  will  be  fully  defined  -functions  of  time  that  prescribe 
the  motion  of  points  on  the  body.  This  permits  the  direct  partitioning  of  the  mass  and  damping 
matrices  in  (4.7). 
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Extracting  only  the  top  equation  for  the  unconstrained  DOF  gives  the  equation  of  motion  with 
boundary  conditions  applied;  that  is. 
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Mq(0  +  Dq(0  +  (q,  t)  =  h(q,  t) 

h(q,  t)  :=  4,  (q,  t)  -  v(0  -  D^v(0 


(4.10) 


Here  the  over  bars  are  used  to  emphasize  that  these  DOF  are  the  free  DOF,  and  the  internal  and 
external  forces  do  depend  on  all  the  DOF.  During  the  linearization  of  no  additional  matrix- 

partition  terms  need  to  be  included  in  h  since  their  contributions  are  already  present  in  f  int  {q,t). 

4.4  Fluid- structure  interactions 
4.4.1  Overview 

There  are  many  varieties  of  partitioned  methods  based  on  a  prediction-correction  model.  The 
methods  differ  in  what  is  used  for  the  prediction,  and  if  the  correction  is  used  repeatedly  or 
staggered  in  time.  The  method  implemented  here,  as  outlined  in  Figure  19,  is  sub-iterated  until 
convergence  of  the  entire  system's  equilibrium  is  achieved.  The  fluid  solver  is  based  on  same 
explicit  fractional  step  method  used,  but  implemented  inside  the  FLASH  framework  (ASC  Flash 
Center,  2012;  Daley  et  ah,  2012).  The  large  scale  high  performance  computing  (HPC) 
framework  is  designed  to  tackle  extremely  large  problem  domains  on  finite-difference  grids.  In 
FLASH,  one  can  use  either  uniform  gridding  techniques  or  adaptive  mesh  refinement  (AMR) 
based  on  the  PARAMESH  library  (MacNeice  et  al.,  2000). 

Each  time  step  is  started  by  predicting  the  states  of  the  structure,  and  computing  the  position, 
velocity,  and  acceleration  fields  of  the  body's  wet  surface.  Then  the  fluid  velocity  field  u  is 


Fractional  Step  DNS 


Figure  19:  Procedure  diagram  for  the  partitioned  FSI  algorithm. 
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computed.  In  the  usual  manner  of  fractional  step  methods,  this  field  is  not  divergence  free.  The 
surface  kinematics  of  the  body  are  then  applied  to  the  fluid  grid  points  around  the  body  in  a 
Lagrangian  fashion  as  implemented  by  Vanella  (  2010).  The  calculation  of  the  pressure  P 
presents  the  most  expensive  step  in  the  calculation.  This  elliptical  problem,  often  referred  to  as 
the  Poisson  problem,  is  discretized  to  become  a  set  a  simultaneous  linear  equations.  The 
efficient  calculation  of  the  pressure  represents  one  of  the  largest  hurtles  to  large  scale  solutions 
(Daley  et  ah,  2012).  Once  the  pressure  gradient  is  computed,  the  corrected  (or  end-of-step) 
velocity  is  calculated  and  stored.  The  velocity  and  pressure  information  are  then  used  to 

compute  the  surface  forces  on  the  body.  A  corrector  procedure  then  computes  an  updated 
estimate  of  the  states  of  the  body.  If  the  states  have  not  changed  within  some  tolerance,  then  the 
velocity  field  of  the  immersed  boundary  conditions  is  recomputed  and  the  cycle  repeats.  Once 
the  convergence  criterion  has  been  satisfied,  time  is  incremented  and  the  outer  loop  begins  again 
with  a  prediction  of  the  structure  using  the  previous  fluid  load. 

The  implementation  issues  of  the  predictor-corrector  method  used  here  are  mostly  surrounding 
the  treatment  of  the  body,  since  the  coding  for  the  fluid  model  was  already  in  place.  This  entails 
both  the  time  integration  of  the  body  as  well  as  the  construction  of  the  forcing  terms  as  boundary 
conditions  on  each  partitioned  field.  Inside  the  FLASH  architecture  there  is  an  entire  unit  of  the 
code  dedicated  to  Lagrangian  particle  tracking  known  as  PARTICLES.  The  previous  uses  of 
these  particles  range  from  physics  simulations,  to  convecting  massless  particles  for  event  tracing. 
For  the  FSI  implementation  considered  here,  they  will  serve  as  the  method  of  communication 
between  the  fluid  domain  and  the  structural  domain. 

4.4.2  Imposition  of  boundary  conditions 

Inside  the  FLASH  code,  the  PARTICLES  unit  is  a  well  apportioned  framework  for  working  with 
Lagrangian  points  distributed  across  the  Eulerian  domain.  The  distribution  of  the  particles  on  the 
HPC  cluster  is  performed  by  FLASH.  The  immersed  boundary  unit  called  ImBound  uses  the 
information  of  each  particle  to  enforce  the  no-slip  condition.  The  use  of  PARTICLES  then  is  to 
cover  the  body's  surface  with  particles  whose  kinematics  is  prescribed  by  the  surface  of  the  body. 
These  particle  points  are  used  for  both  parts  of  the  communication  of  boundary  conditions: 
forcing  the  fluid  and  forcing  the  body. 

In  Figure  20A,  a  small  collection  of  particles  inside  the  fixed  fluid  grid  is  shown.  These  particles 
each  represent  a  small  patch  of  the  surface  area  on  the  surface  of  a  body  (Figure  20B).  It  is  the 
information  of  the  particles  that  connects  the  fluid  and  structural  domains.  Each  patch  of  surface 
must  be  near  the  same  size  as  the  fluid  grid  spacing.  Therefore  the  spacing  of  particles  is 
determined  by  the  fluid  grid  since  for  most  problems  the  fluid  mesh  will  be  much  finer  than  the 
body's  mesh.  This  permits  the  meshing  of  the  structure  to  be  independent  of  the  fluid  grid. 
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(A) 


(B) 


Figure  20:  A  representative  region  of  fluid  particles  in  the  fluid  domain,  and  these  same  particles 
in  the  structural  domain,  (a)  The  fluid  domain  particles  are  used  to  represent  the  kinematics  of  the 
body's  surface,  (b)  Likewise,  the  particles  are  interpreted  hy  the  body  as  patches  of  constant 

applied  traction. 

A  designer  of  the  body  would  need  only  to  be  concerned  with  sufficient  spatial  resolution  for  the 
body's  deformation.  This  is  in  contrast  to  previous  immersed  boundary  implementations.  In  the 
work  of  Vanella,  (2010),  a  rigid  wing  with  the  planform  of  a  Musca  domesica  was  defined  using 
381,662  triangles.  If  each  triangle  was  mapped  to  the  face  of  a  tetrahedron  finite  element,  the 
number  of  DOF  for  a  relatively  simple  body  would  be  staggering.  The  method  implemented 
here  using  particles  avoids  this  complication  by  grouping  particles  through  surface  elements. 
Now  a  subset  of  the  particles  is  indexed  to  a  surface  element,  and  this  mapping  is  structured  to 
allow  for  memory  efficiency  and  calculation  speed  (Fitzgerald,  2013). 

The  mapping  of  each  particle  to  an  individual  patch  of  surface  on  the  body  is  shown  in  Figure 
21.  Here  n  is  the  outward  facing  unit  vector  at  center  of  the  patch,  and  {tj,  tj}  is  a  pair  of 

vectors  tangent  to  the  surface.  The  surface  normal  is  computed  by  the  cross-product  of  two 
independent  vectors  on  the  surface  of  the  body.  To  ensure  that  this  calculation  results  in  outward 
facing  normals,  a  check  is  performed  during  mesh  pre-processing.  The  distributed  force  is 

defined  in  the  global  frame.  The  area  of  each  surface  patch  is  .  The  calculations  to  determine 

the  kinematics  of  the  deformed  surface  are  performed  at  the  center  point  of  the  patch.  Since  the 
spacing  of  particles  is  the  same  as  the  fluid  grid,  then  it  is  assumed  that  fj-  is  constant  on  each 

patch. 
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Figure  21;  Discrete  surface  element  of  the  body.  The  triad  is  a  local  set  of  unit  vectors 

on  the  patch  p ,  with  differential  area  dA  and  surface  traction  f^  . 


The  local  ordering  of  the  patches  on  each  element  is  constructed  by  using  Figure  22.  First  the 
length  of  the  deformed  element  is  calculated,  and  compared  with  the  local  fluid  grid  to  determine 
the  spacing  of  particles  on  the  surface  element.  There  are  and  particles  in  the  ^  -direction 

and  T]  -direction,  respectively.  Each  particle  is  assigned  a  local  index  p  ,  and  this  integer  is  used 
to  uniquely  place  the  patch  in  the  2D  grid  {k,  j)  where  ^  e  [1,  ]  cz  M  and  j  e  [1,  ]  cz  N .  From 

there,  a  unique  mapping  is  defined  to  identify  the  index  p  with  all  the  information  needed  to 
consistently  integrate  the  constant  surface  traction  over  the  patch  and  project  the  forces  to  the 
element’s  degrees  of  freedom.  Also  this  mapping  provides  all  the  necessary  information  to 
reconstruct  the  position,  velocity,  acceleration,  deformed  area,  and  outward  unit-normal  used  to 
impose  the  no-slip  boundary  condition  on  the  fluid  grid. 


V 


Figure  22;  Natural  domain  of  the  surfaee  element,  showing  a  single  particle  patch. 
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4.4.3  The  Generalized-a  method 

The  use  of  highly  resolved  finite  elements  results  in  both  large  storage  requirements  and  also  in 
severe  time  step  requirements  for  explieit  integrators.  Most  finite  element  integrators  use  a 
variant  of  the  seminal  Newmark-|3  method  (Newmark,  1959).  The  Generalized-a  method  (G-a) 
is  a  popular  method  in  the  dynamies  of  linear  problems  dating  baek  to  Chung  and  Hulbert 
(1993).  It  represents  a  unification  of  the  methods  of  Hilber  et  al.  (1977)  and  Wood  et  al.  (1980), 
with  improved  characteristics.  The  G-a  is  second  order  in  time,  implicit,  unconditionally  stable 
for  linear  problems,  and  has  user  selectable  dissipation  of  high  frequencies  .  It  was  shown  to 

be  suitable  for  use  for  nonlinear  problems  in  structural  mechanics  by  Kuhl  and  co-workers  (Kuhl 
and  Crisfield,  1999;  Kuhl  and  Ramm,  1996).  A  major  consequences  of  using  G-a  on  nonlinear 
structures  is  the  loss  of  unconditional  stability.  A  detailed  analysis  of  the  properties  of  the 
method  for  simple  nonlinear  systems  can  be  found  in  the  studies  of  Baldo  et  al.  (2006;  Bonelli  et 
al.  (002),  and  Erlicher  et  al.  (2002). 

The  method  requires  only  slight  modifications  to  handle  essential  boundary  conditions  that  are 
time  dependent.  The  forcing  term  is  modified  per  h  in  (4.10),  and  it  is  interpolated  with  the 
other  external  loads  and  not  with  the  acceleration  terms.  The  implemented  G-a  in  an  implicit 
that  is  consistently  linearized  to  perform  a  full  Newton-Raphson  continuation  of  the  equilibrium 
of  the  dynamic  equations  of  motions.  The  convergence  criterion  for  the  Newton-Raphson 
scheme  was  chosen  to  be  based  on  the  change  to  the  velocities  of  the  degrees  of  freedom.  This 
can  be  stated  as 


error 


7 


kq[ 


(4.11) 


where  7  and  /3  are  constants  determined  by  the  choice  of  .  The  G-a  method  is  easily 

adapted  for  use  in  the  Structure  predictor  and  Structure  corrector  roles  of  Figure  19  with  some 
careful  tweaking.  There  are  three  main  points  where  discretion  is  required  to  use  the  method 
efficiently  and  robustly: 

1 .  The  choice  of  spectral  radius  . 

2.  The  choice  of  the  tolerance  s  in  the  Newton-Raphson  iterations,  and  how  it  should  be 
different  between  a  predictor  step  and  a  corrector  step  in  the  FSI  scheme. 

3.  The  choice  of  FSI  convergence  criterion. 

The  choice  of  the  spectral  radius  is  a  body  specific  problem.  Since  the  radius  is  relative  to  the 
time  step  a  choice  of  a  smaller  time  step,  say  for  the  CFL  condition  of  the  fluid,  would  result  in 
more  temporal  resolution  in  the  structure.  Therefore  knowing  estimates  of  the  time  step  that 
complies  with  the  CFL,  and  knowing  how  many  modes  of  the  body  are  likely  to  be  of  interest 
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provides  an  estimate  for  the  value  of  .  The  bodies  considered  here,  such  as  flapping  wings  or 

moving  plates,  are  mostly  undergoing  bending  deformation  similar  to  their  first  mode  shape. 
The  CFL  condition  at  low  Reynolds  numbers  has  a  much  smaller  At  requirement  over  the  coarse 
FEM  body.  Since  little  excitation  of  the  higher  modes  are  seen,  a  relatively  small  value  of  is 
suitable  for  most  of  these  types  of  problems. 

Changing  the  value  of  s  between  a  prediction  step  and  a  correction  step  was  found  to  be  very 
beneficial  to  the  FSI  convergence  rate.  In  the  first  method  explored  here,  the  predictor  and 
corrector  are  set  to  only  take  a  single  Newton-Raphson  step  that  resulted  in  FSI  substeps  in  the 
40-60  range.  A  critical  review  of  how  the  Newton-Raphson  iterations  take  the  body  from  time  n 
to  «  + 1  provides  the  basis  for  improving  on  that  scheme.  For  bending  problems,  the  first  step  of 
a  Newton-Raphson  method  moves  the  body  in  the  direction  of  bending;  this  makes  sense  since 
that  is  the  direction  with  the  lowest  stiffness.  The  subsequent  Newton-Raphson  shifts  the  body 
axially.  That  first  step  commonly  overestimates  the  displacement,  and  the  subsequent  iterations 
could  be  seen  as  pulling  the  body  back  to  equilibrium.  When  the  single  Newton-Raphson  step 
was  used,  the  convergence  rate  was  impractical  since  the  fluid  was  reacting  to  a  body  that  was 
not  close  to  equilibrium.  Using  these  observations,  the  Newton-Raphson  exit  criterion  s  can  be 
intelligently  changed. 

•  Predictor:  Set  the  value  of  £•  to  lx  10“^.  This  makes  the  predicted  deformation  field  of 
the  body  relatively  close  to  the  previous  deformation. 

•  Corrector:  Set  the  value  of  s  to  lx  10“^.  This  allows  the  Newton-Raphson  several 
substeps  before  returning  to  the  fluid  for  an  updated  set  of  loads.  The  overprediction  of  a 
single  Newton-Raphson  step  is  avoided,  and  the  fluid  loads  are  very  current.  The 
structural  solver  does  not  waste  time  being  very  exact  since  the  fluid  loads  will  change. 

The  FSI  convergence  criterion  must  be  set  such  that  the  change  of  the  states  during  a  correction 
is  less  than  the  tolerance  fpg; ;  that  is,  errorpg,  <  fpgj  where  the  error  is  defined  as 


errorpgi 


(0)^(«+i) 


(4.12) 


In  light  of  how  s  for  the  corrector  is  chosen,  the  condition  on  FSI  must  be  more  strict,  fpgj  <  e . 

Therefore,  a  value  of  £’pg,  =1x10”*  is  chosen  as  the  default.  Also  this  scheme  implies  that  at 

least  one  correction  step  will  always  occur  which  ensures  the  strong  coupling  of  the  equations. 
During  test  simulations  with  the  values  of  s  as  stated  above,  the  number  of  FSI  substeps 
dropped  from  40-60  down  to  5-8  for  the  same  problem  setup.  This  represents  an  incredible 
speed  up  and  opens  the  possibility  of  using  the  method  in  production  simulations. 
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5  Moving  plate  example 

In  order  to  test  the  FSI  code,  a  simple  and  relatively  small  sized  problem  was  required.  Recently, 
Cleaver  et  al.  (2013a,  2013b)  performed  experiments  on  compliant  plates  in  a  water  tunnel. 
Their  force  measurements  could  provide  FSI  numerical  studies  with  a  simple  validation  case. 
The  problem  setup  here  has  been  simplified  to  make  the  domain  manageable  for  a  small  number 
of  nodes  on  the  HPC.  Here,  only  64-128  processors  are  used  with  wall  times  between  12  to  18 
hours. 

5.1  Formulation 

A  moving  plate  of  length  Z,  width  0.3Z,  and  thickness  0.05Z  is  centered  in  a  3Lx3Lx3L 
quiescent  fluid  domain.  The  plate  is  rotated  15  degrees  along  its  long  axis.  The  density  of  the 
body  is  varied  between/?/ e  {1,  2,  10}.  The  displacements  of  one  of  the  short  edges  is 
restrained  by  the  prescribed  kinematics 

Xj  (t)  =  0 

X2(0  =  0  (5.1) 

Xj  (/)  =  (^1  —  j  sin 

where  the  time  constant  r  =  0.05  is  chosen  as  a  small  number  to  remove  an  impulsive  start  but 
not  affect  the  kinematics  for  long.  The  Poisson's  ratio  was  chosen  to  be  0.3  and  Young's 
modulus  was  chosen  such  that  the  first  natural  frequency  of  the  body  was  co^  /  (X>^=\/3 .  The 

maximum  prescribed  velocity  is  chosen  as  the  reference  speed 

=  max  |x3  (t)|  =  (5-2) 

The  amplitude  of  oscillation  is  set  to  ^3  =  0.15Z  ,  while  the  Reynolds  number  is  constructed  as 

(5.3) 

V 

These  parameters  are  not  the  same  as  those  used  by  Cleaver  et  al.  (2013  a,  2013b),  the  Reynolds 
number  has  been  lowered  from  10,000,  there  is  no  free-stream  velocity,  and  the  domain  is 
smaller.  The  setup  tested  here  is  merely  a  first  step  and  the  use  of  future  resources  would  allow 
for  larger  problems  that  exactly  replicates  the  published  work. 
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5.2  Results 


The  solution  is  computed  for  several  periods  7^  =1k  !  of  forcing.  In  Figure  23,  the  surface 

tractions  on  a  body  are  shown.  From  these  plots  it  is  clearly  evident  that  even  for  Re  =  200 ,  the 
viscous  stresses  are  mostly  quite  low,  except  at  the  edges  of  the  body.  The  pressure  is  acting  as 
one  would  expect  with  a  large  high  pressure  opposing  the  motion  of  the  plate. 

Computations  of  the  center  of  mass  show  how  the  effect  of  mass  density  changes  the  FSI  across 
different  configurations.  In  Figure  24,  the  center  of  mass  is  shown  for  densities 
p !  e  {10,  2,  1}  at  Re  =  200  as  well  as  a  dry  body.  The  displacements  in  Xj  and  X2  appear 
to  change  drastically  as  the  density  ratio  is  lowered,  and  the  FSI  forces  of  the  fluid  begin  to 
dominate  the  body.  This  is  especially  visible  in  the  X2  direction  where  the  displacements  are 
nearly  10  times  those  of  the  dry  (non-FSI)  case.  The  fluid  also  appears  to  be  damping  out  the 
higher  frequencies  of  the  response.  This  is  seen  in  the  X3  direction,  where  only  the  dry  body 
appears  to  have  multiple  frequencies  in  the  response. 


Figure  23;  Example  of  the  surfaee  stresses  acting  on  the  body.  These  are  the  instantaneous  values 
at  1 1 =  2.05  ,  for  Re  =  200 ,  p  /  =  1 .  (A)  Pressure  on  the  body,  with  colors  corresponding 

to  C^l  1.  (B)  Magnitude  of  the  viscous  stress,  with  colors  corresponding  to  /  2  . 
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Figure  24:  Position  of  the  center  of  mass  in  the  deformable  plate  for  Re  =  200 . 

Integrating  the  surface  stresses  across  the  body  provides  a  way  to  compare  the  fluid  contributions 
across  the  density  range.  The  total  force  due  to  viscous  stresses  is  quite  small  for  all  time  in 
Figure  25.  This  correlates  well  to  the  instantaneous  Add  shown  in  Figure  23B  where  the  peaks 
may  have  been  large,  but  they  were  highly  localized.  The  pressure  forces  are  dominating  in  all 
three  directions  by  an  order  of  magnitude. 

The  visualization  of  the  flow  fields  of  Figure  26  were  constructed  by  using  TecPlot.  At  several 
instances  of  time,  a  number  of  isocontours  of  the  Q-criterion  (Hunt  et  al.,  1988)  are  plotted  along 
with  a  X2  -  X3  slice  of  the  Xj  velocity  field.  Depicted  in  the  left  column  are  results  for  Re  =  200 , 
and  in  the  right  column  are  Re  =  1000  at  nearly  the  same  instances  of  time.  The  density  ratio  of 
the  body  is  p/  =  1 ,  which  as  demonstrated  in  the  previous  figures  results  in  the  largest 

deformation  and  highest  loads.  The  Q-isosurfaces  are  not  symmetric  about  the  Xj  -  Xj  plane  since 
the  body  is  rotated  along  the  X2-axis.  Comparing  side-by-side  frames  shows  that  the  flow 
structures  appear  smaller  and  remain  longer  as  the  Reynolds  number  is  increased.  At 
t  /  =  1 .61  for  =  1000  the  Q-isosurface  appears  to  be  rolling  up  on  itself  in  a  hairpin-like 
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Figure  25:  Fluid  forces  computed  at  the  centroid  of  the  deformable  plate  for  Re  =  200 . 

manner  (Bernard,  2011).  Overall,  the  results  from  these  tests  appear  promising.  The  FSl 
algorithm  works  efficiently  enough  to  make  serious  problems  practical  for  calculations.  This  is 
assuming  that  the  necessary  HPC  power  is  available. 
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Re  =  200 


Re  =  1000 


Figure  26:  Demonstration  of  FSI  for  a  flexible  plate  with  p  /  =  1 ,  and  co^  !  co^  =  1  /  3  .  The 

time  has  been  nondimensionalized  with  .  Shown  are  isosurfaces  of  the  Q-criterion  and  a  slice 

of  the  z  velocity  in  the  y-z  plane. 
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6  Concluding  Remarks 

6.1  Future  work 

The  next  steps  are  to  validate  the  implemented  FSI  eode  by  using  the  experiments  of  Cleaver  et 
al.  (2013a,  b)  and  use  the  code  to  simulate  flexible  wings.  The  first  insect  inspired  geometry  to 
be  modeled  is  the  forewing  of  a  Manduca  sexta.  In  Figure  27,  the  researchers  show  the  first  draft 
of  a  3D  FEM  mesh  of  the  planform.  This  shape  was  extracted  from  a  figure  provided  in  the 
work  of  O’Hara  and  Palazotto  (  2012).  Another  useful  reference  for  wing  planforms  is  the  work 
of  Comstock  (1918),  which  contains  many  illustrations. 

The  driving  kinematics  can  be  based  on  the  functions  discussed  by  Berman  and  Wang  (2007). 
Their  angular  description  of  flapping  motion  is  applicable  to  a  broad  range  of  flapping  styles. 
The  hawkmoth  parameters  of  Berman  and  Wang  (2007)  suggest  an  almost  constant  angle  of 
attack  for  the  majority  of  the  hover  stroke  with  relatively  fast  reversals  at  the  ends  of  the  stroke. 
This  kinematics  has  been  implemented  in  the  FLASH  code  and  illustrated  in  Figure  28  by  using 
the  parameters  provided  by  Berman  and  Wang  (2007)  for  the  hawkmoth. 


Figure  27:  Mesh  of  a  Manduca  sexta  inspired  wing.  The  planform  is  based  on  the  results  of 
O’Hara  and  Palazotto  (  2012).  The  restrained  nodes  are  at  the  root  in  magenta.  The  free  nodes  are 
shown  in  blue  (only  the  top  surface  is  shown  for  clarity). 
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Figure  28;  Demonstration  of  the  implemented  kinematies  employing  the  hawkmoth  parameters 

from  Berman  and  Wang  (2007). 

Preliminary  dry  tests  with  =  1  /  3  for  a  homogeneous  body  indicate  that  the  body  is  too 

soft.  A  realistic  correction  is  stiffening  of  the  root  and  leading  edge  of  the  wing.  This  follows 
from  the  spirit  of  the  parameter  distribution  proposed  by  Combes  and  Daniel  (2003a)  as  well  as 
the  mode  shapes  found  in  §4.2. 

The  Reynolds  number  of  a  Manduca  sexta  based  on  the  span  length  and  peak  tip  speed  would  be 
in  the  range  of  25x10^ —29x1 0\  This  increase  in  Reynolds  number  drastically  changes  the 
computing  requirements  for  a  practical  simulation.  If  uniform  grids  are  used  for  example, 
assuming  that  the  cell  size  would  need  to  be  approximately  Ax  =  0.002L ,  then  a  small  domain  of 
4T  X  4T  X  4T  would  have  roughly  2048^  points.  If  the  FLASH  block  size  is  32^ ,  this  results  in 
64^  blocks  to  distribute  on  a  HPC.  For  this  sized  blocks,  it  was  found  that  placing  more  than  30 
blocks  per  processor  results  in  poor  scaling.  Therefore,  216,224  blocks  at  30  blocks  per 
processor  results  in  8,739  processors.  That  many  processors  represent  a  major  use  of  resources. 
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and  this  need  is  beyond  the  clusters  currently  available  at  the  University  of  Maryland.  The  use  of 
AMR  and  lower  Reynolds  number  cases  will  need  to  be  considered  as  the  next  practical  steps. 

6.2  Outcomes 

•  Wing  flexibility:  The  models  that  capture  the  physics  of  flapping  flight  in  2D  were 
developed  and  analyzed.  It  was  shown  that  under  assumptions  moderate  flexibility  is 
beneficial. 

•  Model  fidelity:  A  comparison  study  between  DNS  and  UVLM  in  2D  demonstrated  the 
range  of  usefulness  of  the  computationally  inexpensive  UVLM.  This  comparative  study 
indicates  that  the  UVLM  is  well  suited  to  determine  averaged  quantities,  such  as  those 
that  would  be  used  in  an  objective  function  for  optimization  of  design  variables. 
Subsequent  investigations  around  the  design  point  using  the  DNS  would  still  be  required 
to  ensure  the  desired  results. 

•  Insect  wing  experiments:  A  novel  experiment  to  measure  the  vibration  characteristics  of 
a  living  hawkmoth  was  conducted.  These  findings  correlate  well  to  the  work  of  Sims 
(2010)  and  provide  information  to  support  the  1/3  harmonic  forcing  found  in  nature. 

•  High  fidelity  computation:  Since  the  physics  of  flapping  flight  in  3D  is  less  well-known, 
the  uses  of  low-fidelity  models  are  not  yet  reliable.  This  led  the  researchers  to  build  a 
high-performance  FSI  code  built  into  the  FLASH  framework. 

•  Partitioned  FSI  algorithm:  A  novel  partitioned  FSI  algorithm  using  the  Generalized-a 
method  is  described.  This  method  is  suitable  for  large  scale  systems,  with  bodies  of  tens 
of  thousands  of  degrees  of  freedom.  Special  attention  is  given  to  the  time  step 
requirements,  and  how  this  method  decouples  the  fluid  grid  generation  from  the  body 
mesh  generation.  A  consistent  method  to  project  the  boundary  conditions  between  the 
body  and  the  fluid  is  constructed  using  the  Lagrangian  markers  in  the  PARTICLES  unit 
of  FLASH. 

•  Three-dimensional  studies  with  a  flexible  body:  The  highly  flexible  3D  FSI  code  has 
been  demonstrated  as  operational,  and  algorithmically  tuned.  The  numerical 
demonstrations  of  the  above  methods  show  the  effectiveness  of  the  algorithms  and 
implementations.  A  fully  flexible  plate  is  excited  in  a  flow  over  a  range  of  densities  and 
Reynolds  numbers.  It  was  found  that  the  methods  worked  well,  even  for  equal  densities 
between  the  fluid  and  the  solid.  Large  deformations  were  handled  without  issue  by  only 
12  elements,  and  the  FSI  substeps  were  kept  low  (between  5  to  8). 
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